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We present a unified treatment of Zitterbewegung phenomena for a wide class of systems including 
spintronic, graphene, and superconducting systems. We derive an explicit expression for the time- 
dependence of the position operator of the quasiparticles which can be decomposed into a mean 
part and an oscillatory term. The latter corresponds to the Zitterbewegung. To apply our result 
for different systems one needs to use only vector algebra instead of the more complicated operator 
algebra. 
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The Zitterbewegung (ZB) was first regarded as a rela- 
tivistic effect rooted in the Dirac equation and related to 
a 'trembling' or oscillatory motion of the center of a free 
wave packet [lll^. The ZB is caused by the interference 
between the positive and negative energy states in the 
wave packet; the characteristic frequency of this motion 
is determined by the gap between the two states. It was 
believed that the experimental observation of the effect is 
impossible since one would confine the electron to a scale 
of the Compton wavelength H/ttiqc, where toq is the bare 
mass of the electron [2|. However, the ZB is not a strictly 
relativistic effect: it can appear even for a nonrelativis- 
tic particle moving in a crystal [3| or for quasiparticles 
governed by the Bogoliubov-de Gennes equations in su- 
perconductors [J|. 

Most recently, Schliemann et al. [5|,y| predicted the ZB 
in spintronic systems where the experimental observation 
of the effect is more realistic due to the much smaller fre- 
quency of the oscillatory motion. In these semiconductor 
nanostructures [7| spin-orbit coupling generates an os- 
cillatory motion of the wave packet. The semi-classical 
time evolution of holes was investigated numerically for 
the Luttinger Hamiltonian by Jiang et al. [g. The rela- 
tion between the ZB and the spin transverse force was 
studied by Shen [9| . In a numerical work Lee and Bruder 
observed an oscillatory behavior in the charge and spin 
densities of quantum wires with Rashba and Dresselhaus 
types of spin-orbit coupling [lO|. With a spin-polarized 



electron injected into a waveguide, Nikolic et al. [11 1 
found an oscillatory motion of the wave packet numeri- 
cally, and the ZB pattern was also predicted numerically 
by Brusheim and Xu fj]. Similarly, Zawadzki studied 



the ZB in narrow gap semiconductors [13| , in single- wall 
semiconducting carbon nanotubes [14| and in crystals us- 
ing the nearly- free electron approximation [l5|, which is 
essentially the same as the two-band model in [3|. 

Two-dimensional carbon sheets, known as graphene, 
have been studied theoretically [id Il7| for many decades, 
since their band structure is unique, a gapless Dirac-like 



spectrum [l8|. However, the experimental consequences 
of such a relativistic electron dynamics were observed 
only recently in Hall conductivity measurements 19|, |20| . 
In bilayer graphene a more peculiar behavior of the Hall 
effect was observed experimentally 21| , and explained in 
terms of the chiral Hamiltonian first derived by McChann 
and Fal'ko [23. Both in single and bilayer graphene the 
appearance of the oscillatory motion of the electron re- 
lated to the ZB was pointed out by Katsnelson [23| . Most 
recently, Tworzydlo et al. associated the shot noise with 
the interference of electron-hole pairs at the Dirac point 
in graphene 2J| . As an experimental observation of the 
ZB Trauzettel et al. proposed to measure the photon- 
assisted electron transport in graphene 25 1. 

In this work we present a unified description of the ZB 
in the systems mentioned above. Our approach makes it 
possible to calculate with simple algebra (without using 
operator algebra) the time dependence of the position 
operator of the particle for a wide class of systems. We 
also easily verify the results first obtained by Schliemann 
et al. [g, [a| for spintronic systems. Our result directly 
shows that the ZB is not necessarily a relativistic effect 
but it is related to the coupling between the components 
of the eigenstates of the system. This phenomenon is 
thus the direct consequence of the pseudo-spin degree of 
freedom. 

The time-dependence of the position operator in the 
Heisenberg picture is given by r(i) = e'^*/'' r(0) e^*^*/'*, 
where H is the Hamiltonian of the system. To calculate 
the operator r(i) one can work with the eigenstates of 
H. However, a further insight into the nature of the 
ZB can be gained by solving the equations of motion. 
We start with a quite general form of the Hamiltonian 
that is suitable to describe the systems mentioned in the 
introduction: 



H = e{p)t + n^ S, 



(1) 



where the system is characterized by the one-particle 
energy dispersion £(p) and the effective magnetic field 
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TABLE I: The Hamiltonian of different systems can be expressed as in Eq. ((T|. Here D is the dimension of the system, 
P± = Pi ± ipy, a± = ax i: iuy, and the spin operator is S — | cr, where a = {(Jx,Uy,(Jz) is the set of Pauli matrices. For 
Cooper pairs we assume (for simplicity) that the pair potential A is real and independent of r, and that the energy is measured 
from the Fermi energy Ep- In the last row, ep — h^k^ /{2m), where p — Kk, q is fixed, and Vq is the Fourier transform of 
the periodic potential treated as a perturbation in the crystal. Here iR{-} and 3{-} are the real and imaginary parts of the 
argument. More details of these systems can be found in the references listed in the last column. 



ri(p) coupled to the spin S. Here we assume that e(p) 
and ri(p) are differentiable functions of the momentum 
P = {Px ,PyiPz)- Here T stands for the transpose of a vec- 
tor, while 1 is the unit matrix in spin space, which will be 
omitted hereafter. In the absence of an electrostatic po- 
tential V{r) the momentum p, and, consequently, ri(p), 
are constants of motion. In Table U we listed a few sys- 
tems (together with the effective magnetic field J^(p)) 

I 



that are currently intensely studied in spintronics, and 
in the research of graphene and superconductors. 

It should be emphasized that although in Table U the 
Hamiltonian for all systems is given in terms of the Pauli 
matrices corresponding to a spin S = ^, in our general 
consideration, the spin operator S in Eq. ([1]) can repre- 
sent a quasiparticle with an arbitrary spin S > ^. 



r 



The equations of motion of the position operator r(i) 
and the spin operator S(t) in the Heisenberg picture for 
the Hamiltonian ([1]) read 
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H, S] ^n{p) X S(t), where (2b) 
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Note that Eqs. ([2]) are coupled equations of all three com- 
ponents of S{t) and r(i). However, in the case of two- 
dimensional systems only the x- and y-components of r(i) 
are involved in Eq. (|2ap . 

It is clear from Eq. (j2bp that the spin vector S(i) pre- 
cesses around the vector O. The solution of Eq. (j2bp 
with the initial condition S(0) = Sq can be written as 

S(t) = [n o n + (1 - n o n) cos fii + sin ili nx] So, (3) 
where ft — fin, n is a unit vector, fP = fl ft, and 



non denotes the outer or direct product, i.e., (n o n)^^, = 
riink- Here the operator Sq on the right-hand side is in 
the Schrodinger picture, i.e., it is time independent. One 
can show that the usual commutation relations still hold 
[S,{t),Sjit)]^ihe^jkSk{t). 

Inserting Eqs. ^ into Eq. (Pa|) and solving the differ- 
ential equation one finds 
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with the initial condition r(0) = ro. This is our central 
result. The interpretation of the different terms in ^ 
is as follows. The ZB stems from the oscillatory terms 
(cosine and sine terms). In contrast to the usual dynam- 
ics (first and third terms), two new terms appear in the 
non-oscillatory part: the transverse displacement, which 
is independent of time (second term), and a term that 
corresponds to a particle motion with constant anoma- 
lous velocity (third term). In addition to the oscillatory 
part, these two terms in r(i) are inherent of the ZB. The 
anomalous velocity plays a crucial role in the anomalous 
and spin Hall effects in semiconductors [29|. 

To evaluate the time-dependent position operator r{t) 
within a Gaussian wave packet one can follow, e.g., the 
calculation presented in Refs. [5|, la|. 

We are now in a position to apply our results to the 
systems listed in Table ID Using Eqs. (|2c)) and ((4]), some 
simple algebra yields the same results as given by Eqs. (7) 
and (8) in Ref. [6| for the Rashba-Dresselhaus system. 
Similarly, it is easy to verify the results Eqs. (41) and 
(42) in Ref. [6| for systems of heavy holes in a quantum 
well HO]. 

The current operator in graphene systems splits into 
three terms of which the last one can be associated with 
the ZB phenomenon [23| . Our general approach can also 
be applied to graphene layers to find the time evolution 
of the position operator r(t). We now present explicit re- 
sults for the position operator r(i) from which the trem- 
bling (oscillatory) motion of the electron in graphene 
systems is clearly seen. For single-layer and bilayer 
graphene, the explicit formulas for x(t) and y{t) can again 
be easily obtained using Table |T] and Eqs. ([^ and ^. 
The results for single-layer graphene are 
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and for bilayer graphene 
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Here ai are the Pauli matrices and p^ — p^ + Py- 

Similarly, using the Hamiltonian for Cooper pairs given 
in Table [J the following results are obtained: 
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One can show that these results 



Px+Py+Pz 

agree with those presented in Ref. i4|| . 

Similarly, some simple algebra yields the same results 
as in Refs. [3|,ll5| for nearly free electrons listed in Table |I] 
(except that the off-diagonal elements are swapped in 
the latter reference). For bulk Dresselhaus systems (3rd 
row in Table [I]) the calculation is again straightforward 
but the results are rather cumbersome and not presented 
here. 

Discussion. As mentioned above, the position opera- 
tor r(t) in (|4]) is decomposed into a mean part and an 
oscillatory term. If one derives the position operator 
r(i) directly from r{t) = e'"*/^ r(0) e'^"^/^ then such 
a decomposition can only be obtained using the Foldy- 
Wouthuysen transformation [J, |3lj . In this transforma- 
tion the operator r(t) is calculated in the basis of the 
eigenstates of Hamiltonian ([1]). It can be shown that for 
S = \ (with Pauli matrices) the eigenenergies E±{k) and 
the eigenstates V'±(r) = |x±)e''^'" are given by 
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where and $ are the spherical polar angles of the vector 
fl{h'k) in k-space, and |a| is the magnitude of vector a. 



However, for S > ^ the Foldy-Wouthuysen transforma- 
tion is more complicated. The advantage of our approach 
is that it leads directly to the desired decomposition of 
the position operator r(i). 

For pure Rashba coupling and for single-layer graphene 
the ZB can be interpreted as a consequence of the con- 
servation of the total angular momentum J^ = Lz + Sz, 
where L = r x p is the orbital angular momentum (see 
Rcf. [6]). However, in general, Jz is not a constant of 
motion, ie., [H,J] = Jl x S — p x KS ^ 0. As it can 
be readily seen, this is the case, for example, for Rashba- 
Dresselhaus systems where a 7^ and /3 7^ 0, or for heavy 
holes in a quantum well. 

Finally, it should be mentioned that the ZB is related 
to the non-trivial behavior of the conductivity of single 
and bilayer graphenes 23[ since the velocity operator (Pa|) 
does not commute with the Hamiltonian ([T]). The pecu- 
liar behavior of the spin Hall effect may also be related 
to the ZB '9\. 
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